###################################################################################################
## Payson (2021) When Cities Lobby Replication Code
## Chapter Five: Figures (5.1 - 5.7) and Tables (5.1)
## Created with R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out"
###################################################################################################


## Load pacakges
library(maps); library(ggplot2); library(scales)
library(effects); library(stargazer); library(RColorBrewer); 
library(DataCombine); library(lfe); library(tidyverse)
library(fect); library(gutenbergr); library(tidytext)


## Set working directory
setwd()


## Functions for clean tex tables (need to define to compile tables)
pretty <- function(x) {
  paste(prettyNum(x, big.mark = ","))}

num.cities <- function(x) {
  length(which(getfe(x)=="fips"))}

N <- function(x) {
  x$N}

getmean <- function(x) {
  round(mean(x$response), 1)}



###################################################################################################
## Figures 5.1 - 5.7
###################################################################################################

###################################################################################################
## Figure 5.1 Lobbying and State Transfers Across Cities
###################################################################################################

dta <- read.csv("city_lobby_finance.csv")

(p <- ggplot(subset(dta, state.pp1 > 0 & state.pp1 < 300), aes(x = state.pp1, linetype = factor(lobby))) +
   geom_density() + 
   coord_cartesian(xlim=c(-10, 300)) +
   scale_linetype_manual(values=c("solid", "dashed"), name = "", na.translate = FALSE, 
                         labels = c("Non-Lobbying Cities", "Lobbying Cities")) +
   labs(y = "Density", x = "State Transfers Per Capita")  +
   theme_bw(base_family = "serif", base_size = 12) +
   theme(plot.background = element_blank(),
         legend.position = "top",
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank())) 


###################################################################################################
## Figure 5.2 Lobbying and State Transfers: Difference-in-Differences Predictions
###################################################################################################

dta <- read.csv("city_lobby_finance.csv")

m <- felm(state.pp ~ lobby + log.pop | fips + year, data = dta, keepModel = TRUE)


## Get fixed effects and separate by city and year
fe <- getfe(m)

fips.fe <- fe[fe$fe=="fips",]
year.fe <- fe[fe$fe=="year",]


## Calculate predicted values when lobbying = 0 and lobbying = 1
x <- c(0, 1)

y.fit <- weighted.mean(fips.fe$effect, w = fips.fe$obs) + weighted.mean(year.fe$effect, w = year.fe$obs) +
  coef(m)[2]*mean(m$model$log.pop) + coef(m)[1]*x


sse <- sum((m$response - m$fitted.values)^2)
mse <- sse / (m$N - 2)

mean.se.fit0 <- (1 / m$N + (0 - mean(m$model$lobby))^2 / (sum((m$model$lobby - mean(m$model$lobby))^2)))
mean.se.fit1 <- (1 / m$N + (1 - mean(m$model$lobby))^2 / (sum((m$model$lobby - mean(m$model$lobby))^2)))

se0 <- sqrt(mse * mean.se.fit0)
se1 <- sqrt(mse * mean.se.fit1)


plot.dta <- data.frame(effect = y.fit, lower = c(y.fit[1] - 1.9*se0, y.fit[2] - 1.9*se1),
                       upper = c(y.fit[1] + 1.9*se0, y.fit[2] + 1.9*se1),
                       lobby = c(0, 1))

(p <- ggplot(plot.dta, aes(lobby, effect)) +
    geom_point() +
    scale_y_continuous(limits = c(270, 290), 
                       breaks = c(270, 275, 280, 285, 290)) +
    scale_x_continuous(breaks = c(0, 1), limits = c(-.25, 1.25),
                       labels = c("Non-Lobbying Cities", "Lobbying Cities")) +
    geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1) +
    stat_summary(geom="line") +
    theme_bw(base_family = "serif", base_size = 12) +
    theme(plot.margin = unit(c(1,2,1,1), "cm"),
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white")) +
    labs(y = "State Transfers Per Capita ($)", x = ""))



###################################################################################################
## Figure 5.3 Effect of Lobbying on State Transfers: Parallel Trends
###################################################################################################

dta <- read.csv("city_lobby_finance.csv")

## Identify switchers
dta$switchers <- ifelse(dta$lobby==1 & dta$lobby.lag1==0 & dta$lobby.lag2==0, 1, 0)

## Highest number of switchers in 2013 (88)
table(dta$switchers, dta$year)

## Identify cities that didn't lobby for two years and started lobbying in 2013
dta$treat <- ifelse(dta$year==2013 & dta$lobby==1 & dta$lobby.lag1==0 & dta$lobby.lag2==0, 1, 0)


dta <- dta %>%
  group_by(fips) %>%
  mutate(treated = sum(treat))


## Drop extreme outliers and limit to balanced sample
dta1 <- subset(dta, state.pp < 3000 & total.cog == 10)

tmp <- dta1 %>%
  group_by(year, treated) %>%
  summarize(mean.pp = mean(state.pp, na.rm = T),
            median.pp = median(state.pp, na.rm = T),
            log.state = mean(log.state, na.rm = T))


(p <- ggplot(subset(tmp, year > 2008 & year < 2015), aes(year, mean.pp, linetype = factor(treated))) +
    geom_point() + geom_line() +
    scale_linetype_discrete(name = "", labels = c("Non-Lobbying Cities", "Start Lobbying in 2012")) +
    geom_vline(xintercept = 2012.5, linetype = 2) +
    theme_bw(base_family = "serif", base_size = 12) +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white"),
          legend.position = "top") +
    labs(x = "Years Until Treatment", y = "State Transfers Per Capita ($)"))



###################################################################################################
## Figure 5.4 Effect of Lobbying Expenditures on State Transfers
###################################################################################################

dta <- read.csv("city_lobby_finance.csv")

exp.states <- dta[dta$state=="Alaska" | dta$state == "Arizona" | dta$state == "California" | 
                    dta$state == "Colorado" | dta$state == "Connecticut" | dta$state == "Florida" |
                    dta$state == "Georgia" | dta$state == "Idaho" | dta$state == "Illinois" | 
                    dta$state == "Indiana" | dta$state == "Iowa" | dta$state == "Kentucky" |
                    dta$state == "Maine" | dta$state == "Massachusetts" | dta$state == "Michigan" |
                    dta$state == "Mississippi" | dta$state == "Montana" | dta$state == "Nebraska" |
                    dta$state == "New Jersey" | dta$state == "New York" | dta$state == "North Carolina" |
                    dta$state == "Oregon" | dta$state == "Pennsylvania" | dta$state == "South Carolina" |
                    dta$state == "Tennessee" | dta$state == "Texas" | dta$state == "Virginia" | 
                    dta$state == "Washington" | dta$state == "Wisconsin",]


exp.states$amount <- ifelse(exp.states$amount < 1000, 0, exp.states$amount)

exp.states$log.amount <- log(exp.states$amount + 1)
exp.states$lobby.pp <- exp.states$amount / exp.states$pop


m <- lm(state.pp ~ lobby.pp + log.pop + fips + year,data = subset(exp.states, state.pp < 1000))


plot.dta <- data.frame(lobby.pp = c(0, 3.5, 6.9, 10, 14),
                       lower = summary(effect("lobby.pp", m, confidence.level=.9))$lower,
                       upper = summary(effect("lobby.pp", m, confidence.level=.9))$upper,
                       effect = summary(effect("lobby.pp", m, confidence.level=.9))$effect)


(p <- ggplot(subset(plot.dta, lobby.pp <10.1), aes(lobby.pp, effect)) +
    geom_line() +
    geom_line(aes(x = lobby.pp, y = upper), linetype = "dashed") +
    geom_line(aes(x = lobby.pp, y = lower), linetype = "dashed") +
    theme_bw(base_family = "serif", base_size = 12) +
    scale_x_continuous(breaks = seq(0, 10, 2)) +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white")) +
    labs(x = "Lobbying Expenditures Per Capita ($)", y = "State Transfers Per Capita ($)"))



###################################################################################################
## Figure 5.5 Effect of Lobbying by Median Income
###################################################################################################

dta <- read.csv("city_lobby_finance.csv")

m <- felm(state.pp1 ~ lobby*log.inc + log.pop + log.inc + log.own.source + 
              pct.white + log.home  | fips + year, 
            data = dta)

beta.hat <- coef(m)
cov <- vcov(m)
z0 <- seq(min(dta$log.inc, na.rm = T), max(dta$log.inc, na.rm = T), length.out = 1000)
dy.dx <- beta.hat["lobby"] + beta.hat["lobby:log.inc"]*z0
se.dy.dx <- sqrt(cov["lobby", "lobby"] + z0^2*cov["lobby:log.inc", "lobby:log.inc"] + 
                   2*z0*cov["lobby", "lobby:log.inc"])
upr <- dy.dx + 1.9*se.dy.dx
lwr <- dy.dx - 1.9*se.dy.dx

figdata <- data.frame(x = z0, y = dy.dx, min = lwr, max = upr)

(p <- ggplot(figdata, aes(x = x, y = y)) + 
    geom_line() +
    geom_line(aes(x=x, y = upr), linetype = "dashed") +
    geom_line(aes(x=x, y = lwr), linetype = "dashed") +
    xlim(10.5, 12.25) +
    ylim(-25, 125) +
    labs(y = "Effect of Lobbying on State Transfers Per Capita", 
         x = "City Median Income (Log)") +
    theme_bw(base_family = "serif", base_size = 12)  +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white")))



###################################################################################################
## Figure 5.6 Bill Subjects Lobbied in California by City Income
###################################################################################################

dta <- read.csv("ca_bills.csv")

billtext <- data.frame(text = dta$subject, inc = dta$inc.tercile)

tidy_bill <- function(inc) {
  inc %>% 
    unnest_tokens(word, text) %>% 
    anti_join(stop_words)
}

bills <- tidy_bill(billtext) %>%
  mutate(word = str_extract(word, "[:alpha:]+")) %>% 
  count(inc, word, sort = TRUE)

bills <- bills[!is.na(bills$word) & bills$inc!="avg" &
                 bills$word!="afdc" & bills$word!="abc" &
                 bills$word!="ana" & bills$word!="cancer", ]

comparison_df <- bills %>%
  add_count(inc, wt = n, name = "total_word") %>% 
  mutate(proportion = n / total_word) %>% 
  select(-total_word, -n) %>% 
  pivot_wider(names_from = inc, values_from = proportion, 
              values_fill = list(proportion = 0)) %>%
  pivot_longer(3, names_to = "other", values_to = "proportion")

keep <- c("sick", "pharmaceutical", "stewardship", "disposal", "metering", "wireless", 
          "transmission", "emissions", "solar", "paid", "recreational",
          "economic", "agency", "fees", "waste", "air", "contracts", "infrastructure", "facilities", 
          "resources", "program", "agencies", "land", "formula",
          "transportation", "development", "housing", "water", "budget", "act", "local", "cannabis", 
          "zones", "enterprise", "firearms", "gaming", "flood", "urban", "aid", "youth", "literacy", 
          "gang", "ammunition", "revitalization", "hiv", "underperforming", "accessibility", "agriculture", "fines",
          "sewerage", "discrimination", "unemployment", "impounding", "incarceration")

table(keep %in% comparison_df$word)
keep[!(keep %in% comparison_df$word)]

comparison_df$label <- ifelse(comparison_df$word %in% keep, comparison_df$word, NA)

(p <- comparison_df %>% 
    filter(proportion > 1 / 1e5) %>% 
    ggplot(aes(proportion, `rich`)) +
    geom_abline(color = "gray40", lty = 2) +
    geom_jitter(aes(color = abs(`rich` - proportion)),
                alpha = 0.1, size = 1, width = 0.3, height = 0.3) +
    geom_text(aes(label = label), check_overlap = TRUE, vjust = 1.5) + 
    scale_x_log10(labels = label_percent()) +
    scale_y_log10(labels = label_percent()) + 
    xlab("Low Income Cities") +
    ylab(("High Income Cities")) +
    guides(color = FALSE) +
    theme_bw(base_family = "serif", base_size = 12) +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))



###################################################################################################
## Figure 5.7 Median Income and Number of Bills Lobbied in California
###################################################################################################

dta <- read.csv("ca_bills.csv")

bills <- dta %>%
  group_by(city) %>%
  summarize(total.bills = length(city),
            pop = mean(pop, na.rm = T),
            inc = mean(med.inc, na.rm = T),
            log.pop = log(mean(pop, na.rm = T)),
            log.inc = log(mean(med.inc, na.rm = T)))


bills$avg.bills <- bills$total.bills / 14


m <- lm(avg.bills ~ log.pop + log.inc, data = bills)


effect("log.inc", m, confidence.level=.9)

plot.dta <- data.frame(log.inc = c(9.83, 10.5, 11.1, 11.7, 12.3),
                       lower = c(summary(effect("log.inc", m, confidence.level=.9))$lower[1],
                                 summary(effect("log.inc", m, confidence.level=.9))$lower[2],
                                 summary(effect("log.inc", m, confidence.level=.9))$lower[3],
                                 summary(effect("log.inc", m, confidence.level=.9))$lower[4],
                                 summary(effect("log.inc", m, confidence.level=.9))$lower[5]),
                       upper = c(summary(effect("log.inc", m, confidence.level=.9))$upper[1],
                                 summary(effect("log.inc", m, confidence.level=.9))$upper[2],
                                 summary(effect("log.inc", m, confidence.level=.9))$upper[3],
                                 summary(effect("log.inc", m, confidence.level=.9))$upper[4],
                                 summary(effect("log.inc", m, confidence.level=.9))$upper[5]),
                       effect = c(summary(effect("log.inc", m, confidence.level=.9))$effect[1],
                                  summary(effect("log.inc", m, confidence.level=.9))$effect[2],
                                  summary(effect("log.inc", m, confidence.level=.9))$effect[3],
                                  summary(effect("log.inc", m, confidence.level=.9))$effect[4],
                                  summary(effect("log.inc", m, confidence.level=.9))$effect[5]))

limits <- aes(ymin = lower, ymax = upper)


(p <- ggplot(plot.dta, aes(log.inc, effect)) +
    geom_line(colour = "black") +
    geom_line(aes(x = log.inc, y = upper), linetype = "dashed") +
    geom_line(aes(x = log.inc, y = lower), linetype = "dashed") +
    theme_bw(base_family = "serif", base_size = 13) +
    theme(panel.spacing = unit(2, "lines"),
          plot.margin = unit(c(1,2,1,1), "cm"),
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white")) +
    labs(y = "Bills Lobbied Per Year", x = "Median Income (Log)"))



###################################################################################################
## Table 5.1 City Lobbying by Income
###################################################################################################

dta <- read.csv("city_lobby_finance.csv")

dta$inc.cat <- cut(dta$median.income, include.lowest = TRUE, breaks = quantile(dta$median.income, 
                                                        probs = c(0, .25, .5, .75, 1), na.rm = T),
                   labels = c("Bottom 25%", "25th - 50th Percentile", "50th - 75th Percentile", "Top 25%"))


cbind(Pct.Lobbying = round(prop.table(table(dta$inc.cat, dta$lobby), 1)*100)[,2])



